#ifndef CUDA_KERNELS_H
#define CUDA_KERNELS_H
#define PI 3.14159265359f
//Kernels


// SPH Kernels
__device__ float KernelPoly6( float rsqr, float kernelH ) {
	float w = 315.0/(64.0f*PI*(pow(kernelH, 9)));
	float kernelHSqr = kernelH*kernelH;
	if (rsqr < kernelHSqr) {
		return w*pow((kernelHSqr-rsqr), 3.0f);
	}
	
	return 0.0f;
}

__device__ float KernelPoly6Laplac( float r, float kernelH ){
	float w1 = (pow(kernelH, 4.0f));
	float w2 = 6*pow(kernelH, 2.0f);

	float lengsq = r*r;

	if (lengsq > (kernelH*kernelH)){
		return 0.0f;
	}

	return ( -945*(w1 - w2*lengsq + 5*lengsq*lengsq ) / (32*PI*pow(kernelH, 9.0f)) );	
}

__device__ float KernelViscosityLaplac( float r, float kernelH) {
	float w = 45.0/(PI*pow(kernelH, 6.0f));

	if (r > kernelH){
		return 0.0f;
	}

	return w*(kernelH-r);
}

__device__ float KernelPoly6Grad(  float r, float kernelH  ) {

	float w = (32*PI*pow(kernelH, 9.0f));

	float rsq = r*r;
	float kernelHSqr = kernelH*kernelH;
	if (rsq > kernelHSqr){
		return 0.0f;
	}

	return -945 * r * ((kernelHSqr - rsq)*(kernelHSqr-rsq))/w;
}

__device__ 	float KernelSpikeyGrad( float r, float kernelH ){
	float w1 = -45.0f * (1.0f / (PI*pow(kernelH, 6)));

	if (r > kernelH){
		return 0.0f;
	}

	return w1 * pow((kernelH - r), 2);	
}

__device__ 	float KernelSpikey( float r, float kernelH ) {

	if (r > kernelH){
		return 0.0f;
	}

	return (15.0f / (PI*pow(kernelH, 6))) * pow((kernelH - r), 3);
}



#endif